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Abstract 

A program for calculating the semi-classic transport coefficients is described. It 
is based on a smoothed Fourier interpolation of the bands. From this analytical 
representation we calculate the derivatives necessary for the transport distributions. 
The method is compared to earlier calculations, which in principle should be exact 
within Boltzmann theory, and a very convincing agreement is found. 
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LONG WRITE-UP 



1 Introduction 



Method developments, the existence of user friendly distributed codes and the 
ever increasing computer power are making the calculation of band-structures, 
for even relatively complex materials, more and more straight forward. As sev- 
eral properties can be calculated from the energy bands and their derivatives, 
the usefulness off a generally applicable, easily portable and documented code 
for analysis of the bands should be clear. 

The code presented here relies on a Fourier expansion of the band energies 
where the space group symmetry is maintained by using star functions. The 
idea of the Fourier expansion is to use more star functions than band en- 
ergies, but to constrain the fit so that the extrapolated energies are exactly 
equal to the calculated band-energies and use the additional freedom to mini- 
mize a roughness function and thereby suppress oscillations between the data- 
points. (0; 12 0) Using the analytical representation of the bands it is then a 
reasonable simple procedure to calculate band-structure dependent quantities. 
The method has been tested for several applications based on Boltzmann the- 
ory, including the transport coefficients of intermetallic compounds, (0) high 
T c superconductors(5) m d thermoelectric*. © Furthermore the preset code 
has already been applied to calculate the transport coefficients in a series of 
different clathrate structures^ and a very good agreement was found with ex- 
perimental values. The good agreement was also found for the demanding 
Hall coefficient that depends on the second derivative of the bands. (0; 0) 

Because of the known limitations of Boltzmann theory (fioh the comparison 
with experimental measurements is not the best method for testing the actual 
algorithm for expanding the bands. As the interpolated bands pass exactly 
though the calculated band energies, the precision of the present method is 
mainly limited by possible band crossings where the band derivatives will 
be calculated wrongly. We will therefore test our method by comparing with 



the resent results by Scheidemantel et al.(|lll). Scheidemantel et al. calcu- 



lated transport coefficients of Bi 2 Te 3 (JllJ) by calculating the group velocities 



from the momentum matrix elements. As the momentum matrix elements can 
be calculated directly from the wavef unction. |l2h their method should avoid 



any problems at band crossings. ([ill) As the calculations were documented in 
detail ([ 111) and Bi 2 Te 3 has a complex band-structure that is strongly influenced 
by spin orbit coupling, it constitutes a challenging test-case which we will use 
in the present paper. 
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2 Code implementation 



2.1 Algorithms 



The code relies on a Fourier expansion of the band energies where the space 
group symmetry is maintained by using star functions 

e 4 (k) = £piufifc(k) , S R (k) = -£e* kAR (l) 

R 71 {A} 

where R is a direct lattice vector, {A} are the n point group rotations. The 
idea of the Fourier expansion is to use more star functions than band energies, 
but to constrain the fit so £$ are exactly equal to the band-energies, and use 
the additional freedom to minimize a roughness function. (0; The choice 
of the roughness function, pr, was discussed by Pickett et al.(|3j) who found 
the following expression to be useful for suppressing oscillations between the 
data-points. 

'* = 1 1 - Cl f iicri) +c *(w L ) (2) 

where R m j n is smallest nonzero lattice vector. C\ and Ci are parameters, 
but our and earlier^ experience found that the results are quite insensitive 
to their actual value and we have therefore fixed them to C\ = Ci — 3/4. 
To ensure that £j pass exactly through the calculated the band energies at 
the same time as the roughness function is minimized, the algorithm needs 
sufficient freedom. This means that the number of planewaves must be larger 
than the number of band energies. The number of planewaves to the number 
of band-energies is controlled by the input parameter LPFAC, Table 1, and the 
program prints a warning if the fit is poor (subroutine KCOMP). 

The expansion coefficients are given 

c m = <H kjv) " Er ^° ^ CRi6lk ' R R = ° (3) 

where Ar are calculated by solving 

Ae<(k) = e<(k) - £i (k N ) = J2 #kk'A R (4) 

kVkiv 



where 



^kk' = 2^ I 5 ) 

R^O PR 
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The time consuming steps in the Fourier expansion are obviously the solution 
of Eq. (4) and the setup of the i?kk' matrix, Eq. (5). The matrix setup can 
be identified as a multiplication of two k x R matrices and both Eq. (4) and 



Eq. (5) can thus be handled efficiently by BLAS calls. (113t Il4j) The calculation 



of the expansion parameters, cr,^, are carried out in the subroutine FITE4. 
3 Test problems 

3.1 Boltzmann theory: The semi-classic equations 
Boltzmann 

theory fi3 El G3) 

is a useful tool for gaining insight into the trans- 
port properties of real materials. In the presence of an electric and magnetic 
field and a thermal gradient the electric current, j, can be written in terms of 
the conductivity tensors 

ji = o-ijEj + o-ijkEjBk + fijVjT H (6) 

In terms of the group velocity 



and the inverse mass tensor 



the conductivity tensors can be obtained 

o- aP (i, k) = e 2 T i:k v a (i, k)^, k) (9) 
while a a [3 y is most elegantly written using the Levi-Civita symbol, e^t.|l6Hl7l) 

a a/3y (i, k) = e 3 rl k e^ uv v a (i, k)v v (i, k)M^ (10) 

The notation used in Eqs. (9-10) gives directly the symmetry of the conduc- 
tivity tensors. F.inst. in an orthorhombic symmetry o a p is diagonal with all 
three components independent and <j Q/ g 7 has three independent components 
and vanishes unless a, (3 and 7 are all different. 

The relaxation time, r, in principle is dependant on both the band index 
and the k vector direction. However detailed studies of the direction de- 
pendence of r have shown that, to a good approximation, r is direction 
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independent (18) and that even in the superconducting cuprates, that have 
substantially anisotropic conduction and cell-axes, the r is almost isotropic, (jl) 
In the present we will use the simplest approximation for the relaxation time, 
namely to keep it constant, which is the most often used in praxis. 



Similar to the density of states energy projected conductivity tensors can be 
defined using the conductivity tensors, Eqs. (9-10) 



where N is the number of k-points sampled. Similarly cr aj a 7 (e) can be defined. 
With the expansion of the bands, Eq. (1), the necessary derivatives, Eq. (7), are 
straightforwardly calculated as Fourier sums which can be efficiently evaluated 
using fast Fourier transforms (FFTs). Evaluation of the density of states and 
transport distributions thus requires a total of 10 FFTs for each band in the 
general case. The calculation of the transport distributions is carried out in 
the subroutine DOS and are output to the files: case.transdos, case.sigxx, 
case.sigxxx. 

The transport tensors, Eq. (6), can then be calculated from the conductivity 
distributions 



dMT-e) 



eTQ 
1 



a a p(e)(e - p) 



de 

de 
9U(T;e) 



cr a p(e)(£ - /if 



de . 
df,(T;e 



a af 3 7 {T; n) 



n 



dfJT;e) 



de 



de 
de 



de 



de 



(12) 
(13) 
(14) 
(15) 



where k° is the electronic part of the thermal conductivity. The Seebeck and 
Hall coefficients can then easily be calculated 



Si, = E^T)- 1 = (a- 1 ) 



Hjk = — J — ■ = ( rr 



Rink — — — (c 1 ) Q jO" a/ 3fc((T 1 ) i/ 3 



■appl rjappl 
Ji D k 



(16) 
(17) 



Under the assumption that the relaxation time r is direction independent, 
both the Seebeck and the Hall coefficients are independent of r. The integrals 
Eqs. (12)-(15) are performed in the subroutine FERMIINTEGRALS. 
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(b) 




Fig. 1. Transport coefficients as a function of chemical potential: (a) Seebeck coeffi- 
cient (b) Power factor with respect to scattering time S 2 a/r. One obtains the power 
factor in the usual units of //W/(cm K 2 ) by multiplying by r in units of 10~ 14 s. 



T 300 K ..... / 


V X *- 




df/de 
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0.2 



0.3 



Fig. 2. Integrand factor df/de in Eqs. (12)-(15) at T = 300 K. Arbitrary y units. 
3.1.1 Test case: Bi 2 Te3 

The calculation was carried out using the WIEN code(J_9) with the same 
computational parameters as in ref. (|llh. The calculated transport coefficients 
were found to be converged using a non-shifted mesh with 56000 k points (4960 
in the IBZ). The original fc-mesh was interpolated onto a mesh four times as 
dense. 

The calculated transport coefficients are given in Figure 1. Figure la,b show 
the Seebeck coefficient and the power factor with respect to scattering time. 
Both these curves can be compared to the earlier work (jTlh and an excellent 
agreement is found, both with respect to shape and absolute values. 

Fig. 1 demonstrates that potential problems at band crossings have negligible 
influence on the calculated transport coefficients for Bi2Te3 at 300 K. It should 
off-course be underlined that this is just one example and systems could exist 
where the present method should fail. However, band crossings only happen on 
symmetry lines (symmetry planes in hexagonal systems), so, while the problem 
exists, as long as the k-sampling is dense enough to keep the error localized 
at the crossing it will have little effect on global quantities like transport. 
Furthermore, at T = 300 K, the df/de factor, Fig. 2 is quite broad and has 
5 % of its maximum value ate — /i = 0.11 eV. The transport coefficients are 
thus a sum over several Fermi surfaces and any problems will be smeared out. 
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Fig. 3. Electronic thermal conductivity as a function of chemical potential. Full line: 
calculated from Eqs. (14) and (18). Dashed line: Calculated from Wiedemann- Franz 
law, Eq. (19) 

As a further illustration and test of the method we have calculated the elec- 
tronic thermal conductivity at zero electric current. This can be defined as 

4 = k° - Tv ioi {a- l )p a v Pj (18) 

or aproximated through the electronic conductivity via Wiedemann- Franz law, 
which for degenerate charge carriers is given as 

< - t(t)^ t < 19 > 

Figure lc shows the electronic thermal conductivity calculated with two dif- 
ferent methods, Eqs. (18) and (19). As expected the two lines are in very 
good agreement. The second term in Eq. (18), which is directly related to the 
power factor (Fig. 1), is obviously insignificant far from the band-gap, where 
the Seebeck coefficient is small and the conductivity large. Close to the band 
gap it is a significant correction, as illustrated in the small insert in Fig. 3. 



3.1.2 Test case: CoSb 3 

Calculations on CoSb3 were carried out using the Engel-Vosko GGA.(20) The 
unit cell was Im3 with a = 9.0385 A. The Sb atom is placed at the g Wyck- 
off position with (0,0.33537,0.15788). The plane-wave (PW) cut-off was de- 
fined by min(i?A4"T) m ax(fc n )=5.5 corresponding to approximately 588 PW. The 
Brillouin-zone (BZ) was sampled on a shifted tetrahedral mesh with 300 k 
points (17 in the IBZ) for the self consistent calculation. For the transport 
calculations a non-shifted mesh with 24000 k points (1030 in the IBZ) was 
used. The necessary derivatives were then calculated on a FFT grid five times 
as dense. 



The band structure of the skutterudite CoSb3 has been subject of some dis- 
cussion and is sensitive to the lattice parameter and the exchange correlation 
function. (21; 2^; 2^) It is generally agreed that it has parabolic bands close to 
the Fermi level, which we also find, Fig. 4. In this region of parabolic bands 
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[e/u.c] [e/u.c] 



Fig. 4. Band-structure of C0SD3 together with the inverse Hall coefficient (1/Rh and 
the calculated number of carriers n (density of states) and the difference between 
1/Rh and n. 



the Hall coefficient should be inversely proportional to the number of carriers. 
The calculation of the Hall coefficient depends on the second derivatives of the 
bands, Eq. (10), and therefore serves as a demanding test of the precision of 
the method and CoSb3 has therefore been chosen as a test example. Figure 4 
illustrates how 1/Rh and the doping are almost equivalent in the region of 
parabolic bands, while they differ when the region of flat bands. 



4 Input parameters 



Table 1 gives the input parameters used, bandstyle gives the format of the 
band-structure input. The present version of the code is interfaced to the band- 
structure code WIEN2k(|19i). but can easily be interfaced to any other band- 
structure codes, de defines how fine the mesh for the conductivity distribution 
should be, Eq. (11). ecut defines the range of bands used around ef erm in the 
integrals, Eqs. (12)-(15). setgap and gapsize can be used to apply a scissors 
operator to force a certain band-gap. lpfac defines how much denser the 
interpolated mesh should be and thereby the R-cutoff in Eq. 1. The programs 
outputs the conductivity tensors on a grid of T and /1, Eqs. (12)-(15), defined 
by ef cut, tmax and deltat. All output of the program is in Si-units 
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Table 1 

Input variables. 



bandstyle 


Format of band structure to be input 


i debug 


controls level of output 


ef erm 


Fermilevel (in Ry) 


deltae 


de (in Ry) 


ecut 


cut-off energy around Fermilevel 


setgap 


logic switch for band gap manipulation 


gapsize 


new band-gap (in Ry) 


lpf ac 


# of times the interpolated mesh should be denser than the calculated 


ef cut 


range of \i in which the integrations should be performed 


tmax 


max temperature at which the integrations should be performed 


deltat 


temperature step 



5 Porting the code 



The present version of the code is interfaced to the band-structure code 
WIEN2k. (jl9l ) However, as the method uses only the crystal structure and the 
eigen-energies on a mesh as data the code is very easy to interface to other 
band-structure codes. The necessary crystal structure and band-structure in- 
formation is contained in the MODULE band-structure and the user should 
therefore only supply a subroutine that sets up the module. 



6 Conclusion 



We have implemented and tested a method for obtaining an analytical repre- 
sentation of the band-structure. We have applied it to the calculation of trans- 
port coefficients. The method has been compared with an earlier calculation (jlll 
which in principle should be exact within Boltzmann theory, and we found a 
very convincing agreement. 

It should be pointed out that the present method also has several advantages. 
First of all, when an analytical expression of the bands is found, they can be in- 
terpolated onto a finer /c-mesh. Secondly, because only the energies are needed 
the code is easily portable to any band-structure code and furthermore does 
not require the storage of potentially large wavefunction files on disk. Finally, 
second derivatives necessary for the Hall coefficient are straightforwardly cal- 
culated which is not straightforward from the wavefunction itself. Even higher 
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Table 2 

The variables in MODULE bandstructure. 



aac_dir (3 , 3) 


T"* T7* AT / o\ 

REAL (8) 


Direction cosines lor the conventional direct 


aacj:ec(3,3) 




and reciprocal unit cell 


p2c_dir(3,3) 


REAL (8) 


Conversion matrix tor the primitive to conventional 


p2c_rec(3,3) 




unit cell conversion 


nsym 


INTEGER 


number of symmetry operators 


symop(3,3,48) 


T"» T — 1 AT / Ci \ 

REAL (8) 


symmetry operators with respect 






to the direct primitive lattice 


nband 


1JN 1EGER 


number ot bands 


nkpt 


INTEGER 


number of k-points in the IBZ 


xkpoint ( : , : ) 


REAL (8) 


k-points in basis of primitive reciprocal lattice vectors. 






Should be allocated as xkpoint (3, nkpt) 


bandenergy ( : , : ) 


REAL (8) 


eigen-energies in Ry. 






Should be allocated as bandenergy (nband , nkpt ) 



derivatives, necessary e.g. for the calculation of magneto-resistance, could eas- 
ily be calculated with the present method, but it remains to be seen whether 
the accuracy is high enough. 
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